Spatial Workshop at Kelantan State Health Department

Spatial Workshop at Kelantan State Health Department

Trend, Incidence, Choropleth Map and Spatial Autocorrelation

1.0 Prepare environment

library(sf)

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.4.1

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(janitor)

Warning: package 'janitor' was built under R version 4.4.1

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(gtsummary)
library(DT)
library(stringr)
library(readxl)
library(broom)
library(tmap)

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')

library(mapview)

Warning: package 'mapview' was built under R version 4.4.1

library(lubridate)
library(gganimate)

Warning: package 'gganimate' was built under R version 4.4.1

library(spdep)

Warning: package 'spdep' was built under R version 4.4.1

Loading required package: spData

Warning: package 'spData' was built under R version 4.4.1

To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

2.0 Load Data

2.1 Base Map / Polygon Data

kel <- st_read("kelantan.shp")

Reading layer `kelantan' from data source 
  `C:\Users\ACER\Desktop\spatial workshop\kelantan.shp' using driver `ESRI Shapefile'
Simple feature collection with 66 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
Projected CRS: Kertau (RSO) / RSO Malaya (m)

2.2 Population Data

kel_mukim <- read_xlsx("pop_kel.xlsx")
kel_mukim %>% datatable()

2.3 Merging population data to polygon

kel_map <- merge(kel,kel_mukim,by.x="MUKIM", by.y="MUKIM", all.x=T, sort=F)

3.0 Plot Kelantan Map

st_geometry(kel_map)

Geometry set for 66 features 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
Projected CRS: Kertau (RSO) / RSO Malaya (m)
First 5 geometries:

POLYGON ((485501.8 669698.8, 485717.3 669694.6,...

POLYGON ((487716.5 665649.5, 487615.4 665445.1,...

POLYGON ((482744.8 660223.4, 482823.6 660137.8,...

POLYGON ((486936.9 677358.5, 486990.5 677333.7,...

POLYGON ((490841.7 668783.4, 490906.1 668691, 4...

3.1a Kelantan according to Mukim

plot(kel_map[,1]) # Mukim

3.1b Interactive Map

mapview(kel_map[,1])

3.2 Map the population

kel_map %>% ggplot() + 
  geom_sf(aes(fill = JUM_JANTIN)) + 
  ggtitle('Population of Kelantan') + 
  theme_bw()

mapView(kel_map, zcol = "JUM_JANTIN", layer.name = "Population", popup = kel_map$MUKIM)

4.0 Leptospirosis Data

lepto <- read_xlsx("leptospirosis.xlsx") %>% clean_names()
glimpse(lepto)

Rows: 1,106
Columns: 21
$ diagnosis              <chr> "Leptospirosis", "Leptospirosis", "Leptospirosi…
$ tahun_daftar           <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ epid_daftar            <dbl> 6, 9, 11, 12, 18, 18, 19, 17, 21, 23, 34, 34, 3…
$ age                    <dbl> 30, 23, 39, 43, 31, 34, 48, 55, 48, 39, 24, 45,…
$ alamat_semasa_kejadian <chr> "FELDA ARING", "LADANG U&I CIKU 8", "SYARIKAT S…
$ poskod                 <dbl> 18300, 18300, 18300, 18300, 18300, 18300, 18300…
$ latitude_rso           <dbl> 478031, 459494, 441802, 488502, 496760, 441782,…
$ longitude_rso          <dbl> 548141, 564966, 547551, 547701, 539072, 547566,…
$ latitude_wgs           <dbl> 4.944824, 5.034372, 4.897434, 4.945540, 4.93648…
$ longitude_wgs          <dbl> 102.3491, 102.1477, 101.9605, 102.3498, 102.454…
$ notifikasi_status      <chr> "Daftar Kes", "Daftar Kes", "Daftar Kes", "Daft…
$ race                   <chr> "Foreigner", "Foreigner", "Foreigner", "Foreign…
$ kewarganegaraan        <chr> "Bukan Warganegara", "Bukan Warganegara", "Buka…
$ gender                 <chr> "Male", "Male", "Male", "Male", "Male", "Male",…
$ nationality            <chr> "INDONESIA", "INDONESIA", "INDONESIA", "INDONES…
$ klasifikasi_kes        <chr> "Sporadic", "Sporadic", "Sporadic", "Sporadic",…
$ cara_pengesanan_kes    <chr> "Pasif", "Pasif", "Pasif", "Pasif", "Pasif", "P…
$ jenis_rawatan          <chr> "Wad Perubatan", "Jabatan Kecemasan & Kemalanga…
$ daerah                 <chr> "GUA MUSANG", "GUA MUSANG", "GUA MUSANG", "GUA …
$ mukim                  <chr> "CHIKU", "CHIKU", "GALAS", "CHIKU", "CHIKU", "B…
$ negeri                 <chr> "KELANTAN", "KELANTAN", "KELANTAN", "KELANTAN",…

4.1 Convert Leptospirosis data to spatial data

Use st_as_sf ( ) function to convert line listing data that contained Lat/Long attributes to spatial object

lepto <- st_as_sf(lepto, 
                    coords = c("longitude_wgs", "latitude_wgs"), 
                    crs = 4326)
lepto %>% datatable()

Now, Leptospirosis data will be completed with geometry data for each row/case.

Confirm CRS is wgs84

st_crs(lepto)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

4.2 Convert shapefile to RSO to match with Kelantan Map

lepto_2 <- st_transform(lepto, 3168)
lepto_2 %>% datatable()

4.3 Plot map to see outlier

ggplot() +
  geom_sf(data = lepto_2) +
  ggtitle("Map of Leptospirosis") +
  theme_bw()

4.4 Select point only in Kelantan (all_kel)

all_kel <- lepto_2 %>% 
  mutate(within_kel_map = lengths(st_within(lepto_2, kel_map)))
all_kel2 <- all_kel %>% 
  filter(within_kel_map == 1)

5.0 Plot the Cases

5.1 Overall Plot

overall_plot <- ggplot() +
  geom_sf(data = kel) +   #base map
  geom_sf(data = all_kel2, color = "black", size = 0.5) +  #cases in spatial data
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +  
  theme(plot.title = element_text(size = 12),  strip.text = element_text(size = 12)) # cosmetic

overall_plot 

5.1a Plot yearly Leptospirosis Cases

overall_plot + 
  facet_wrap(~tahun_daftar) + #to plot according to year
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5)) +
  theme(plot.title = element_text(size = 12),  strip.text = element_text(size = 12)) +
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022")

5.1b Plot Leptospirosis Cases according to districts

overall_plot + 
  facet_wrap(~daerah) + #to plot according to daerah
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5)) +
  theme(plot.title = element_text(size = 12),  strip.text = element_text(size = 12)) +
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022")

5.2a Animation of yearly Leptopirosis cases

overall_plot_animate <- ggplot() +
  geom_sf(data = kel) +   # base map
  geom_sf(data = all_kel2, color = "red") +  # cases in spatial data
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +  
  theme(plot.title = element_text(size = 12), 
        strip.text = element_text(size = 12)) + 
  transition_time(tahun_daftar) +  # animate over the year'tahun_daftar')
  labs(subtitle = "Year: {frame_time}")  # display the current year in the animation

# To animate the plot
animate(overall_plot_animate, nframes = 7, fps = 1)

5.2b Animation of Leptopirosis cases according to Districts

overall_plot_animate_district <- ggplot() +
  geom_sf(data = kel) +   # base map
  geom_sf(data = all_kel2, color = "red") +  # cases in spatial data
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +  
  theme(plot.title = element_text(size = 12), 
        strip.text = element_text(size = 12)) + 
  transition_states(daerah, transition_length = 2, state_length = 1) +  # animate over the daerah (districts)
  labs(subtitle = "Daerah: {closest_state}")  # display the current daerah in the animation

# To animate the plot
animate(overall_plot_animate_district, nframes = 10, fps = 1)

6.0 Trend of cases

6.1 Yearly Trend

lepto_year <- all_kel2 %>% 
  group_by(tahun_daftar) %>% 
  count() %>% 
  print(n = 7)

Simple feature collection with 7 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 383744.8 ymin: 510209.5 xmax: 503623.9 ymax: 689075
Projected CRS: Kertau (RSO) / RSO Malaya (m)
# A tibble: 7 × 3
  tahun_daftar     n                                                    geometry
*        <dbl> <int>                                            <MULTIPOINT [m]>
1         2016   361 ((388781.6 517176.5), (392434.5 553262.1), (396233.1 52195…
2         2017   158 ((404206.5 576231.2), (420431.3 559879.2), (421662.5 62454…
3         2018   177 ((385416.8 516312), (393507.1 546418.7), (397873.6 523539.…
4         2019    97 ((386545.8 527808.1), (389174 515964.7), (393525.1 546380.…
5         2020    42 ((383744.8 510209.5), (388995.8 515828.9), (395908 522469.…
6         2021    38 ((435521.9 602877.8), (440599.7 536450.7), (442034.2 54136…
7         2022   233 ((388876.4 516777.7), (395312.6 545653.5), (395346.8 54561…

ggplot(lepto_year, aes(x=tahun_daftar, y=n)) +
  geom_line(size = 1) +
  labs(x="Year", y="Number of Cases ", title="Yearly Leptospirosis Cases (2016-2022)") +
 theme_bw()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

6.2 Weekly Trend

lep_week <- all_kel2 %>% 
  group_by(tahun_daftar, epid_daftar) %>% 
  count() %>% 
  print(n = 52)

Simple feature collection with 277 features and 3 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 383744.8 ymin: 510209.5 xmax: 503623.9 ymax: 689075
Projected CRS: Kertau (RSO) / RSO Malaya (m)
# A tibble: 277 × 4
   tahun_daftar epid_daftar     n                                       geometry
 *        <dbl>       <dbl> <int>                                 <GEOMETRY [m]>
 1         2016           1     4 MULTIPOINT ((461636.4 665564.8), (465307.3 68…
 2         2016           2     7 MULTIPOINT ((453598 628501), (454513.6 670146…
 3         2016           3    17 MULTIPOINT ((433218.4 611282.1), (462483.8 58…
 4         2016           4     4 MULTIPOINT ((438269.6 541530.5), (465339.2 67…
 5         2016           5     2 MULTIPOINT ((466310 654978.5), (472532 677163…
 6         2016           6     5 MULTIPOINT ((401048.5 567447.6), (439631 5355…
 7         2016           7     4 MULTIPOINT ((415953 631506), (473923.2 686361…
 8         2016           8     3 MULTIPOINT ((466508 629432), (469531.9 637935…
 9         2016           9     9 MULTIPOINT ((441738.2 539696), (447146 574459…
10         2016          10     6 MULTIPOINT ((442691.3 665699.6), (454915.2 63…
11         2016          11     8 MULTIPOINT ((396233.1 521954.8), (441093.2 54…
12         2016          12     8 MULTIPOINT ((424238 542459), (441945.8 540279…
13         2016          13     6 MULTIPOINT ((455928.7 666323.5), (456702.6 63…
14         2016          14     8 MULTIPOINT ((429155.9 540578.1), (443335.8 66…
15         2016          15    12 MULTIPOINT ((439820.8 536510.8), (446482.8 57…
16         2016          16     8 MULTIPOINT ((441945.8 540279.3), (445116.7 54…
17         2016          17    11 MULTIPOINT ((404692.7 554047.8), (414301.4 54…
18         2016          18     9 MULTIPOINT ((443196.8 567564.1), (458105.7 66…
19         2016          19     8 MULTIPOINT ((419261.9 577001), (424118 542451…
20         2016          20     5 MULTIPOINT ((442748.7 563911.6), (456532.6 66…
21         2016          21     3 MULTIPOINT ((455737.7 688343.3), (470265.4 67…
22         2016          22     2 MULTIPOINT ((439120 529842), (458709.4 683913…
23         2016          23     5 MULTIPOINT ((388781.6 517176.5), (446161.2 57…
24         2016          24     5 MULTIPOINT ((432957.1 622502.3), (458705.2 68…
25         2016          25     6 MULTIPOINT ((421769 542219.1), (448502.2 6344…
26         2016          26     7 MULTIPOINT ((442258 541173), (443794.5 665502…
27         2016          27     1                        POINT (443264.6 665142)
28         2016          28     7 MULTIPOINT ((456754 665421), (458734.9 672812…
29         2016          29     1                        POINT (448033 665595.3)
30         2016          30     6 MULTIPOINT ((400993.7 567424.1), (417903.8 58…
31         2016          31     3 MULTIPOINT ((446482.8 572005.5), (467067 5513…
32         2016          33     7 MULTIPOINT ((392434.5 553262.1), (433100.4 63…
33         2016          34     9 MULTIPOINT ((428179 540245), (449419.6 626755…
34         2016          35     2 MULTIPOINT ((434758.1 645424.5), (460491.2 67…
35         2016          36    11 MULTIPOINT ((424579.9 542693.7), (438141.2 65…
36         2016          37     4 MULTIPOINT ((444834.9 543522.2), (446208.9 65…
37         2016          38    17 MULTIPOINT ((440295 649643.1), (440880.4 6600…
38         2016          39     9 MULTIPOINT ((434787.3 633913.3), (448198 6636…
39         2016          40    17 MULTIPOINT ((441202 640528), (442618 665057),…
40         2016          41     9 MULTIPOINT ((428322.1 629253), (428552.5 6294…
41         2016          42     6 MULTIPOINT ((468359.8 648377.3), (470140.1 63…
42         2016          43     9 MULTIPOINT ((417437 635971), (429544.5 627791…
43         2016          44     6 MULTIPOINT ((441175.2 538192), (450345.3 6390…
44         2016          45    13 MULTIPOINT ((416860 632551), (439752 535440),…
45         2016          46     5 MULTIPOINT ((451643.3 672253.4), (456532.6 66…
46         2016          47     8 MULTIPOINT ((417363 636079), (425457.4 631230…
47         2016          48    12 MULTIPOINT ((446476 548113), (450151.8 662635…
48         2016          49    14 MULTIPOINT ((441869 540303), (442473.7 637940…
49         2016          50     2 MULTIPOINT ((440034.7 638088.2), (458013.8 59…
50         2016          51     7 MULTIPOINT ((446345.7 657884.3), (461212.7 68…
51         2016          52     4 MULTIPOINT ((461779 678961.1), (480949.1 5951…
52         2017           1     1                          POINT (458842 570731)
# ℹ 225 more rows

ggplot(lep_week, aes(x=epid_daftar, y=n, group=tahun_daftar, color = as.factor(tahun_daftar))) +
  geom_line(size = 0.8) +
  labs(x="Epid Week", y="Number of Cases", title="Leptospirosis according to Epid Week and Year", color="Year") +
  scale_color_brewer(palette="Set2") + 
  scale_y_continuous(labels=function(x) format(x, digits=1, nsmall=0)) +
 theme_bw()

ggplot(lep_week %>% filter(tahun_daftar == "2022"), 
       aes(x=epid_daftar, y=n)) +
  geom_line(size = 0.8, color = "blue") +
  labs(x="Epid Week", y="Number of Cases", title= " Leptopsirosis in 2022") +
  scale_color_brewer(palette="Set2") + 
  scale_y_continuous(labels=function(x) format(x, digits=1, nsmall=0)) +
  theme_bw()

6.3 Animation of weekly Leptopirosis cases in 2022

ggplot(lep_week %>% filter(tahun_daftar == "2022"), 
       aes(x = epid_daftar, y = n)) +
  geom_line(size = 0.8, color = "blue") +
  labs(x = "Epid Week", y = "Number of Cases", title = "Leptospirosis in 2022") +
  scale_y_continuous(labels = function(x) format(x, digits = 1, nsmall = 0)) +
  theme_bw() +
  transition_reveal(epid_daftar)  # Animate based on 'epid_daftar' (epidemiological week)

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

7.0 Aggregated Data

For this exercise, we focus our analysis on Leptospirosis cases reported in 2016 only.

lepto_16 <- lepto_2 %>% filter(tahun_daftar == "2016")

7.1 Joint point to polygon

#lepto density per mukim
lepto_muk <- st_join(lepto_16, kel_map, 
                      join = st_within)
glimpse(lepto_muk)

Rows: 361
Columns: 34
$ diagnosis              <chr> "Leptospirosis", "Leptospirosis", "Leptospirosi…
$ tahun_daftar           <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ epid_daftar            <dbl> 6, 9, 11, 12, 18, 18, 19, 17, 21, 23, 34, 34, 3…
$ age                    <dbl> 30, 23, 39, 43, 31, 34, 48, 55, 48, 39, 24, 45,…
$ alamat_semasa_kejadian <chr> "FELDA ARING", "LADANG U&I CIKU 8", "SYARIKAT S…
$ poskod                 <dbl> 18300, 18300, 18300, 18300, 18300, 18300, 18300…
$ latitude_rso           <dbl> 478031, 459494, 441802, 488502, 496760, 441782,…
$ longitude_rso          <dbl> 548141, 564966, 547551, 547701, 539072, 547566,…
$ notifikasi_status      <chr> "Daftar Kes", "Daftar Kes", "Daftar Kes", "Daft…
$ race                   <chr> "Foreigner", "Foreigner", "Foreigner", "Foreign…
$ kewarganegaraan        <chr> "Bukan Warganegara", "Bukan Warganegara", "Buka…
$ gender                 <chr> "Male", "Male", "Male", "Male", "Male", "Male",…
$ nationality            <chr> "INDONESIA", "INDONESIA", "INDONESIA", "INDONES…
$ klasifikasi_kes        <chr> "Sporadic", "Sporadic", "Sporadic", "Sporadic",…
$ cara_pengesanan_kes    <chr> "Pasif", "Pasif", "Pasif", "Pasif", "Pasif", "P…
$ jenis_rawatan          <chr> "Wad Perubatan", "Jabatan Kecemasan & Kemalanga…
$ daerah                 <chr> "GUA MUSANG", "GUA MUSANG", "GUA MUSANG", "GUA …
$ mukim                  <chr> "CHIKU", "CHIKU", "GALAS", "CHIKU", "CHIKU", "B…
$ negeri                 <chr> "KELANTAN", "KELANTAN", "KELANTAN", "KELANTAN",…
$ geometry               <POINT [m]> POINT (484205.4 546893.9), POINT (461889 …
$ MUKIM                  <chr> "CHIKU", "CHIKU", "GALAS", "CHIKU", "CHIKU", "B…
$ NEGERI                 <chr> "KELANTAN", "KELANTAN", "KELANTAN", "KELANTAN",…
$ DAERAH                 <chr> "GUA MUSANG", "GUA MUSANG", "GUA MUSANG", "GUA …
$ LELAKI                 <int> 14709, 14709, 19644, 14709, 14709, 12006, 12006…
$ PEREMPUAN              <dbl> 11384, 11384, 17311, 11384, 11384, 11135, 11135…
$ JUM_JANTIN             <dbl> 26093, 26093, 36955, 26093, 26093, 23141, 23141…
$ `2016`                 <dbl> 32000, 32000, 47800, 32000, 32000, 31700, 31700…
$ `2017`                 <dbl> 32000, 32000, 47800, 32000, 32000, 31700, 31700…
$ `2018`                 <dbl> 32500, 32500, 48800, 32500, 32500, 32700, 32700…
$ `2019`                 <dbl> 32900, 32900, 49800, 32900, 32900, 33600, 33600…
$ `2020`                 <dbl> 33300, 33300, 50800, 33300, 33300, 34600, 34600…
$ `2021`                 <dbl> 33300, 33300, 50800, 33300, 33300, 34600, 34600…
$ `2022`                 <dbl> 33300, 33300, 50800, 33300, 33300, 34600, 34600…
$ average                <dbl> 4571.429, 4571.429, 6828.571, 4571.429, 4571.42…

7.2 Count all leptospirosis cases for each mukim in 2016

count_lep_mukim_yr <- lepto_muk %>% 
  count(DAERAH, MUKIM, tahun_daftar,average) %>% 
  ungroup()
count_lep_mukim_yr %>% datatable()

7.3 Calculate incidence of leptospirosis per 1000 population for mukim in 2016

count_lep_muk_y_1000 <- count_lep_mukim_yr %>% 
  mutate(incidence_lep = (n/average)*1000)
count_lep_muk_y_1000 %>% datatable()

8.0 Plot incidence map Leptospirosis reported in 2016

8.1 Joining the incidence data to basemap of Kelantan

kelmap_with_lep <- st_join(kel_map, count_lep_muk_y_1000)

8.2 Map the count of cases in 2016 for each Mukim

kelmap_with_lep %>% ggplot() + 
  geom_sf(aes(fill = n)) +
  ggtitle('Count of Leptospirosis 2016') + 
  theme_bw() 

# Create a custom popup with multiple fields
popup_info <- paste0(
  "<strong>MUKIM: </strong>", kelmap_with_lep$MUKIM.x, "<br>",
  "<strong>DAERAH: </strong>", kelmap_with_lep$DAERAH.x, "<br>",
  "<strong>Leptospirosis Cases: </strong>", kelmap_with_lep$n
)

# Display the map with the custom popup
mapView(kelmap_with_lep, zcol = "n", layer.name = "Count of Leptospirosis cases 2016", popup = popup_info)

8.3 Map the incidence of cases in 2016 for each Mukim

kelmap_with_lep %>% ggplot() + 
  geom_sf(aes(fill = incidence_lep)) +
  ggtitle('Count of Leptospirosis 2016') + 
  theme_bw() 

# Create a custom popup with multiple fields
popup_info_1<- paste0(
  "<strong>MUKIM: </strong>", kelmap_with_lep$MUKIM.x, "<br>",
  "<strong>DAERAH: </strong>", kelmap_with_lep$DAERAH.x, "<br>",
  "<strong>Leptospirosis Incidence : </strong>", kelmap_with_lep$incidence_lep
) 

mapview(kelmap_with_lep,zcol = "incidence_lep", layer.name = "Incidence of  Leptospirosis cases 2016", popup = popup_info_1)

9.0 Spatial Autocorrelation

Spatial autocorrelation is used to describe the extent to which a variable is correlated with itself through space. This concept is closely related to Tobler’s First Law of Geography, which states that “everything is related to everything else, but near things are more related than distant things” (Tobler 1970). Positive spatial autocorrelation occurs when observations with similar values are closer together (i.e., clustered). Negative spatial autocorrelation occurs when observations with dissimilar values are closer together (i.e., dispersed). 

9.1 Set Neighbouring Polygon for autocorrelation analysis (using QUEEN approach)

Read data Leptospirosis for 2016

linelist <- read_xlsx("leptospirosis.xlsx")
LEP_kel16 <- filter(linelist, `Tahun Daftar` == "2016")

#Moran's I for leptospirosis 2016
polyLEP_kel16 <- merge(kel_map, LEP_kel16, by = c("MUKIM"))
count_lep16 <- polyLEP_kel16 %>% 
  count(DAERAH, MUKIM, average) %>% 
  ungroup()
nb_lep16 <- poly2nb(count_lep16, queen = TRUE) #set neigbouring queen
nb_lep16[[1]]

[1] 17 40

9.2 Assign Weight for Mukim that contribute to Case in observed mukim

lw_lep16 <- nb2listw(nb_lep16, style = "W" , zero.policy = TRUE) #assign weight
lep16_lag <- lag.listw(lw_lep16, count_lep16$n) #create lag function

9.3 Global Moran’s I

moran.test(count_lep16$n, lw_lep16)

    Moran I test under randomisation

data:  count_lep16$n  
weights: lw_lep16    

Moran I statistic standard deviate = 4.8685, p-value = 5.623e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.394961860      -0.018518519       0.007213068 

9.4 Moran Plot

moran.plot(count_lep16$n, lw_lep16, main = "Moran Plot for Leptospirosis 2016")

9.5 Local Moran’s I

We have seen that the Global Moran’s II provides an index to assess the spatial autocorrelation for the whole study region. There is often interest in providing a local measure of similarity between each area’s value and those of nearby areas. Local Indicators of Spatial Association (LISA) (Anselin 1995) are designed to provide an indication of the extent of significant spatial clustering of similar values around each observation. 

Here, we use the localmoran() function to compute the Local Moran’s II for the housing prices data. We set alternative = "greater" which corresponds to testing H0H0: no or negative spatial autocorrelation vs. H1H1: positive spatial autocorrelation

Local_Moran_lep16 <- localmoran(count_lep16$n,lw_lep16)
LEP16_poly <- cbind(count_lep16, Local_Moran_lep16)

Then, we identify the clusters of each type by using the information provided by the Moran’s II scatterplot obtained with the moran.plot() function 

mp <- moran.plot(as.vector(scale(count_lep16$n)), lw_lep16)

We create the variable quadrant denoting the type of cluster for each of the areas using the quadrant corresponding to its value and its spatially lagged value, and the p-value. Specifically, areas with quadrant equal to 1, 2, 3, and 4 correspond to clusters of type high-high, low-low, high-low, and low-high, respectively. Areas with quadrant equal to 5 are non-significant.

m_lep16 <- count_lep16$n - mean(count_lep16$n)# Center the variable of interest around its mean

m_local_lep16 <- Local_Moran_lep16[,1]- mean(Local_Moran_lep16[,1]) #Center the local Moran around the mean


signif <- 0.05 #significant threshold

# builds a data quadrant
quadrant_lep16 <- vector(mode="numeric",length = nrow(Local_Moran_lep16))
quadrant_lep16 [m_lep16>0 & m_local_lep16>0] <- 4
quadrant_lep16 [m_lep16<0 & m_local_lep16<0] <- 1
quadrant_lep16 [m_lep16<0 & m_local_lep16>0] <- 2
quadrant_lep16 [m_lep16>0 & m_local_lep16<0] <- 3
quadrant_lep16 [Local_Moran_lep16[,5]] <- 0

9.6 Plot the LISA Map

LISA_16 <- cbind(LEP16_poly, quadrant_lep16)

LISA_Map_16 <- tm_shape(kel_map) + #base map
  tm_polygons(col = "#f7f7f7")+ 
  tm_style("natural") + 
tm_shape(LISA_16) +
  tm_polygons(col = "quadrant_lep16", breaks = c(0,1,2,3,4,5), palette= 
                              c("white", "blue","skyblue","lightpink","red"),labels = c("Non Significant","Low-Low","Low-High","High-Low", "High-High"), title="LISA 2016") +
  tm_layout(legend.outside = TRUE)

LISA_Map_16

# Define a color palette for the LISA quadrants
colors <- c("blue", "skyblue", "lightpink", "red")

# Plot the map using mapview
LISA_Map_16_mapview <- mapview(LISA_16, 
          zcol = "quadrant_lep16",  # column for color mapping
          col.regions = colors,     # apply the color palette
          legend = TRUE,            # display the legend
          layer.name = "LISA 2016", # name of the layer in the viewer
          at = c(0, 1, 2, 3, 4), # breaks for quadrants
          labels = c("Non Significant", "Low-Low", "Low-High", "High-Low", "High-High"))

Warning: Found less unique colors (4) than unique zcol values (5)! 
Interpolating color vector to match number of zcol values.

# Display the map
LISA_Map_16_mapview